setwd(out_path)
if(!dir.exists("PartySupport")){
	stop("Run model estimation first")
}else{
	setwd("PartySupport")
	if(!dir.exists("post")){
		dir.create("post")
	}else{
		warning("Overwriting previous run, if any")
	}
}
setwd("post")

###############
#Minmax sequence
###############
minmax_seq<-function(numvec){
	.out<-seq(from=min(numvec),
		to=max(numvec),
		by=(diff(range(numvec)/10))
	)
	return(.out)
}

############
#Calculate expected value, with confidence intervals, for POLR
############
polr_expval<-function(rep.full.op,sim.data,aic_clean,id){
	#Calculate expected value, with confidence intervals, for POLR
	
	#Coef matrix
	coefMCraw<-coef(rep.full.op)
	# coefMC2rm<-grep("factor(GovtControl)",names(coefMCraw),value=T,fixed=T)
	# coefMC2rm<-coefMC2rm[!grepl("Dem",coefMC2rm,fixed=T)]
	
	# coefMC<-coefMCraw[setdiff(names(coefMCraw),coefMC2rm)]
	coefMC<-coefMCraw
	coefMC.mat<-matrix(coefMC,nrow=length(coefMC),dimnames=list(names(coefMC)))
	
	#X matrix
	X.mat<-as.matrix(sim.data[,c(names(coefMC)[1:(length(coefMC))])])
	# X.mat[,"GovtControl"]<-1
	class(X.mat)<-"numeric"
	
	#Calculate linear predictor
	stopifnot(dim(X.mat)[2]==dim(coefMC.mat)[1])
	Xb<-X.mat%*%coefMC.mat
	
	#Calculate cutpoints
	cutpointsMC.full<-rep.full.op$zeta
	cutpointsMC<-matrix(cutpointsMC.full,nrow=1)
	
	#Get confidence intervals
	SEs<-as.matrix(sqrt(diag(X.mat%*%vcov(rep.full.op)[names(coefMC),names(coefMC)]%*%t(X.mat)))*1.96)
	pred.df<-NULL
	levs<-as.numeric(rep.full.op$lev)
	Pprev<-0
	
	for(i in levs[1:(length(levs))]){
		ind<-which(levs==i)
		
		#Add in the predictions
		stopifnot(ind<=(length(cutpointsMC)+1))
		if(ind<=length(cutpointsMC)){
			lcur<-cutpointsMC[,ind]-Xb
			lcur.lo<-lcur-SEs
			lcur.hi<-lcur+SEs
			pLeqCur<-pnorm(lcur)-Pprev
			pLeqCurlo<-pnorm(lcur.lo)-Pprev
			pLeqCurhi<-pnorm(lcur.hi)-Pprev
			Pprev<-pnorm(lcur)
		}else if(ind==(length(cutpointsMC)+1)){
			lcur<-qnorm(1-Pprev)
			lcur.lo<-lcur-SEs
			lcur.hi<-lcur+SEs
			pLeqCur<-pnorm(lcur)
			pLeqCurlo<-pnorm(lcur.lo)
			pLeqCurhi<-pnorm(lcur.hi)
		}
		
		#Create dataframe
		cur.df<-sim.data
		cur.df$Rsup<-i
		cur.df$fit<-pLeqCur
		cur.df$lo<-pLeqCurlo
		cur.df$hi<-pLeqCurhi
		pred.df<-rbind(pred.df,cur.df)
	}
	pred.df$id<-id
	return(pred.df)
	
}

###############
#Predicting Republican support from Affluent vs. middle class support
###############

#Fake Data for Graphs (Affluent Opinion)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = minmax_seq(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_aff<-polr_expval(rep.full.op,sim.data,aic_clean,"Affluent Opinion")
ohat_aff$scaled_op<-ohat_aff$pred90_sw
ohat_aff[,c("pred50_sw","pred90_sw")]<-NULL

#Fake Data for Graphs (Middle Class Opinion)
sim.data<-data.frame(b.scaled.intgrp =mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp),
                     pred50_sw = minmax_seq(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_mc<-polr_expval(rep.full.op,sim.data,aic_clean,"Middle Class Opinion")
ohat_mc$scaled_op<-ohat_mc$pred50_sw
ohat_mc[,c("pred50_sw","pred90_sw")]<-NULL

#Now Let's melt the ordered df
ohat.full<-rbind(ohat_mc, ohat_aff)
ohat.full$Rsup<-factor(ohat.full$Rsup,levels=seq(-2,2,1),
					labels=c("Strongly against","Against","Neutral","For","Strongly for"))

#Middle Class Graph
post.rep.sup.mc <- ggplot(ohat.full[ohat.full$id=="Middle Class Opinion",])+
							geom_line(aes(x=scaled_op,y=fit))+
							geom_ribbon(aes(x=scaled_op,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Middle Class Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Republican attitude across middle class opinion",
								subtitle="By agreement level")

#Affluent Graph
post.rep.sup.aff <- ggplot(ohat.full[ohat.full$id=="Affluent Opinion",])+
							geom_line(aes(x=scaled_op,y=fit))+
							geom_ribbon(aes(x=scaled_op,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Affluent Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Republican attitude across affluent opinion",
								subtitle="By agreement level")

###############
#Predicting Democratic support from Affluent vs. middle class support
###############

#Fake Data for Graphs (Affluent Opinion)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = minmax_seq(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_aff<-polr_expval(dem.full.op,sim.data,aic_clean,"Affluent Opinion")
ohat_aff$scaled_op<-ohat_aff$pred90_sw
ohat_aff[,c("pred50_sw","pred90_sw")]<-NULL

#Fake Data for Graphs (Middle Class Opinion)
sim.data<-data.frame(b.scaled.intgrp =mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp),
                     pred50_sw = minmax_seq(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_mc<-polr_expval(dem.full.op,sim.data,aic_clean,"Middle Class Opinion")
ohat_mc$scaled_op<-ohat_mc$pred50_sw
ohat_mc[,c("pred50_sw","pred90_sw")]<-NULL

#Now Let's melt the ordered df
ohat.full<-rbind(ohat_mc, ohat_aff)
ohat.full$Rsup<-factor(ohat.full$Rsup,levels=seq(-2,2,1),
					labels=c("Strongly against","Against","Neutral","For","Strongly for"))

#Middle Class Graph
post.dem.sup.mc <- ggplot(ohat.full[ohat.full$id=="Middle Class Opinion",])+
							geom_line(aes(x=scaled_op,y=fit))+
							geom_ribbon(aes(x=scaled_op,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Middle Class Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Democratic attitude across middle class opinion",
								subtitle="By agreement level")

#Affluent Graph
post.dem.sup.aff <- ggplot(ohat.full[ohat.full$id=="Affluent Opinion",])+
							geom_line(aes(x=scaled_op,y=fit))+
							geom_ribbon(aes(x=scaled_op,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Affluent Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Democratic attitude across affluent opinion",
								subtitle="By agreement level")


demsupGrid<-arrangeGrob(post.dem.sup.aff,post.dem.sup.mc,nrow=2)
repsupGrid<-arrangeGrob(post.rep.sup.aff,post.rep.sup.mc,nrow=2)
dev.off()

ggsave("post_dem_class.tiff",demsupGrid)
ggsave("post_rep_class.tiff",repsupGrid)

###########################
#Predicting Democratic support from Business vs. adv group preferences
###########################

#Fake Data for Graphs (Affluent Opinion)
sim.data<-data.frame(b.scaled.intgrp = minmax_seq(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_aff<-polr_expval(dem.full.op,sim.data,aic_clean,"Bus IntGrp Alignment")
ohat_aff$scaled_ig<-ohat_aff$b.scaled.intgrp
ohat_aff[,c("b.scaled.intgrp","m.scaled.intgrp")]<-NULL

#Fake Data for Graphs (Middle Class Opinion)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = minmax_seq(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_mc<-polr_expval(dem.full.op,sim.data,aic_clean,"Adv IntGrp Alignment")
ohat_mc$scaled_ig<-ohat_mc$m.scaled.intgrp
ohat_mc[,c("b.scaled.intgrp","m.scaled.intgrp")]<-NULL

#Now Let's melt the ordered df
ohat.full<-rbind(ohat_mc, ohat_aff)
ohat.full$Rsup<-factor(ohat.full$Rsup,levels=seq(-2,2,1),
					labels=c("Strongly against","Against","Neutral","For","Strongly for"))

#Middle Class Graph
post.dem.sup.bus <- ggplot(ohat.full[ohat.full$id=="Bus IntGrp Alignment",])+
							geom_line(aes(x=scaled_ig,y=fit))+
							geom_ribbon(aes(x=scaled_ig,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Bus IntGrp Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Democratic attitude across business opinion",
								subtitle="By agreement level")

#Affluent Graph
post.dem.sup.adv <- ggplot(ohat.full[ohat.full$id=="Adv IntGrp Alignment",])+
							geom_line(aes(x=scaled_ig,y=fit))+
							geom_ribbon(aes(x=scaled_ig,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Adv IntGrp Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Democratic attitude across advocacy group opinion",
								subtitle="By agreement level")

###########################
#Predicting Republican support from Business vs. adv group preferences
###########################

#Fake Data for Graphs (Affluent Opinion)
sim.data<-data.frame(b.scaled.intgrp = minmax_seq(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_aff<-polr_expval(rep.full.op,sim.data,aic_clean,"Bus IntGrp Alignment")
ohat_aff$scaled_ig<-ohat_aff$b.scaled.intgrp
ohat_aff[,c("b.scaled.intgrp","m.scaled.intgrp")]<-NULL

#Fake Data for Graphs (Middle Class Opinion)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = minmax_seq(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 0,
					 econ.policies = 0)
ohat_mc<-polr_expval(rep.full.op,sim.data,aic_clean,"Adv IntGrp Alignment")
ohat_mc$scaled_ig<-ohat_mc$m.scaled.intgrp
ohat_mc[,c("b.scaled.intgrp","m.scaled.intgrp")]<-NULL

#Now Let's melt the ordered df
ohat.full<-rbind(ohat_mc, ohat_aff)
ohat.full$Rsup<-factor(ohat.full$Rsup,levels=seq(-2,2,1),
					labels=c("Strongly against","Against","Neutral","For","Strongly for"))

#Middle Class Graph
post.rep.sup.bus <- ggplot(ohat.full[ohat.full$id=="Bus IntGrp Alignment",])+
							geom_line(aes(x=scaled_ig,y=fit))+
							geom_ribbon(aes(x=scaled_ig,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Bus IntGrp Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Republican attitude across business opinion",
								subtitle="By agreement level")

#Affluent Graph
post.rep.sup.adv <- ggplot(ohat.full[ohat.full$id=="Adv IntGrp Alignment",])+
							geom_line(aes(x=scaled_ig,y=fit))+
							geom_ribbon(aes(x=scaled_ig,ymin=lo,ymax=hi),alpha=0.15)+
							facet_wrap(~Rsup)+
							theme_bw()+xlab("Adv IntGrp Opinion")+
							ylab("Probability")+
							ggtitle("Ordered probit, Republican attitude across advocacy group opinion",
								subtitle="By agreement level")

demsupGrpGrid<-arrangeGrob(post.dem.sup.bus,post.dem.sup.adv,nrow=2)
repsupGrpGrid<-arrangeGrob(post.rep.sup.bus,post.rep.sup.adv,nrow=2)
dev.off()

ggsave("post_dem_grp.tiff",demsupGrpGrid)
ggsave("post_rep_grp.tiff",repsupGrpGrid)

setwd(script_path)